/*
 *  iaf_delta_noleak.cpp
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

/* iaf_delta_noleak is a neuron where the potential jumps on each spike arrival. */

#include "iaf_delta_noleak.h"
// Includes from libnestutil:
#include "numerics.h"

// Includes from nestkernel:
#include "exceptions.h"
#include "kernel_manager.h"
#include "universal_data_logger_impl.h"

// Includes from sli:
#include "dict.h"
#include "dictutils.h"
#include "doubledatum.h"
#include "integerdatum.h"

#include <limits>
namespace nest
{

/* ---------------------------------------------------------------- 
   * Recordables map
   * ---------------------------------------------------------------- */

RecordablesMap<iaf_delta_noleak> iaf_delta_noleak::recordablesMap_;

// Override the create() method with one call to RecordablesMap::insert_()
// for each quantity to be recorded.
template <>
void RecordablesMap<iaf_delta_noleak>::create()
{
  // use standard names whereever you can for consistency!
  insert_(names::V_m, &iaf_delta_noleak::get_V_m_);
}

/* ---------------------------------------------------------------- 
 * Default constructors defining default parameters and state
 * ---------------------------------------------------------------- */

nest::iaf_delta_noleak::Parameters_::Parameters_()
    : tau_m_(10.0) // ms
      ,
      c_m_(250.0) // pF
      ,
      t_ref_(2.0) // ms
      ,
      E_L_(-70.0) // mV
      ,
      I_e_(0.0) // pA
      ,
      V_th_(-55.0 - E_L_) // mV, rel to E_L_
      ,
      V_min_(-std::numeric_limits<double>::max()) // relative E_L_-55.0-E_L_
      ,
      V_reset_(-70.0 - E_L_) // mV, rel to E_L_
      ,
      with_refr_input_(false), subtraction_(false)
{
}

nest::iaf_delta_noleak::State_::State_()
    : y0_(0.0),
      y3_(0.0),
      r_(0),
      refr_spikes_buffer_(0.0)
{
}

/* ---------------------------------------------------------------- 
 * Parameter and state extractions and manipulation functions
 * ---------------------------------------------------------------- */

void nest::iaf_delta_noleak::Parameters_::get(DictionaryDatum &d) const
{
  def<double>(d, names::E_L, E_L_); // Resting potential
  def<double>(d, names::I_e, I_e_);
  def<double>(d, names::V_th, V_th_ + E_L_); // threshold value
  def<double>(d, names::V_reset, V_reset_ + E_L_);
  def<double>(d, names::V_min, V_min_ + E_L_);
  def<double>(d, names::C_m, c_m_);
  def<double>(d, names::tau_m, tau_m_);
  def<double>(d, names::t_ref, t_ref_);
  def<bool>(d, "refractory_input", with_refr_input_);
  def<bool>(d, "subtraction", subtraction_);
}

double nest::iaf_delta_noleak::Parameters_::set(const DictionaryDatum &d)
{
  // if U0_ is changed, we need to adjust all variables defined relative to U0_
  const double ELold = E_L_;
  updateValue<double>(d, names::E_L, E_L_);
  const double delta_EL = E_L_ - ELold;

  if (updateValue<double>(d, names::V_reset, V_reset_))
    V_reset_ -= E_L_;
  else
    V_reset_ -= delta_EL;

  if (updateValue<double>(d, names::V_th, V_th_))
    V_th_ -= E_L_;
  else
    V_th_ -= delta_EL;

  if (updateValue<double>(d, names::V_min, V_min_))
    V_min_ -= E_L_;
  else
    V_min_ -= delta_EL;

  updateValue<double>(d, names::I_e, I_e_);
  updateValue<double>(d, names::C_m, c_m_);
  updateValue<double>(d, names::tau_m, tau_m_);
  updateValue<double>(d, names::t_ref, t_ref_);

  if (V_reset_ >= V_th_)
    throw BadProperty("Reset potential must be smaller than threshold.");

  if (c_m_ <= 0)
    throw BadProperty("Capacitance must be >0.");

  if (t_ref_ < 0)
    throw BadProperty("Refractory time must not be negative.");

  if (tau_m_ <= 0)
    throw BadProperty("Membrane time constant must be > 0.");

  updateValue<bool>(d, "refractory_input", with_refr_input_);
  updateValue<bool>(d, "subtraction", subtraction_);

  return delta_EL;
}

void nest::iaf_delta_noleak::State_::get(DictionaryDatum &d, const Parameters_ &p) const
{
  def<double>(d, names::V_m, y3_ + p.E_L_); // Membrane potential
}

void nest::iaf_delta_noleak::State_::set(const DictionaryDatum &d, const Parameters_ &p, double delta_EL)
{
  if (updateValue<double>(d, names::V_m, y3_))
    y3_ -= p.E_L_;
  else
    y3_ -= delta_EL;
}

nest::iaf_delta_noleak::Buffers_::Buffers_(iaf_delta_noleak &n)
    : logger_(n)
{
}

nest::iaf_delta_noleak::Buffers_::Buffers_(const Buffers_ &, iaf_delta_noleak &n)
    : logger_(n)
{
}

/* ---------------------------------------------------------------- 
 * Default and copy constructor for node
 * ---------------------------------------------------------------- */

nest::iaf_delta_noleak::iaf_delta_noleak()
    : Archiving_Node(),
      P_(),
      S_(),
      B_(*this)
{
  recordablesMap_.create();
}

nest::iaf_delta_noleak::iaf_delta_noleak(const iaf_delta_noleak &n)
    : Archiving_Node(n),
      P_(n.P_),
      S_(n.S_),
      B_(n.B_, *this)
{
}

/* ---------------------------------------------------------------- 
 * Node initialization functions
 * ---------------------------------------------------------------- */

void nest::iaf_delta_noleak::init_state_(const Node &proto)
{
  const iaf_delta_noleak &pr = downcast<iaf_delta_noleak>(proto);
  S_ = pr.S_;
}

void nest::iaf_delta_noleak::init_buffers_()
{
  B_.spikes_.clear();   // includes resize
  B_.currents_.clear(); // includes resize
  B_.logger_.reset();   // includes resize
  Archiving_Node::clear_history();
}

void nest::iaf_delta_noleak::calibrate()
{
  B_.logger_.init();

  const double h = Time::get_resolution().get_ms();

  V_.P33_ = std::exp(-h / P_.tau_m_);
  V_.P30_ = 1 / P_.c_m_ * (1 - V_.P33_) * P_.tau_m_;

  // TauR specifies the length of the absolute refractory period as
  // a double in ms. The grid based iaf_psp_delta can only handle refractory
  // periods that are integer multiples of the computation step size (h).
  // To ensure consistency with the overall simulation scheme such conversion
  // should be carried out via objects of class nest::Time. The conversion
  // requires 2 steps:
  //     1. A time object r is constructed defining  representation of
  //        TauR in tics. This representation is then converted to computation time
  //        steps again by a strategy defined by class nest::Time.
  //     2. The refractory time in units of steps is read out get_steps(), a member
  //        function of class nest::Time.
  //
  // The definition of the refractory period of the iaf_delta_noleak is consistent
  // the one of iaf_neuron_ps.
  //
  // Choosing a TauR that is not an integer multiple of the computation time
  // step h will leed to accurate (up to the resolution h) and self-consistent
  // results. However, a neuron model capable of operating with real valued spike
  // time may exhibit a different effective refractory time.
  //

  V_.RefractoryCounts_ = Time(Time::ms(P_.t_ref_)).get_steps();
  assert(V_.RefractoryCounts_ >= 0); // since t_ref_ >= 0, this can only fail in error
}

/* ---------------------------------------------------------------- 
 * Update and spike handling functions
 */

void nest::iaf_delta_noleak::update(Time const &origin,
                                    const long from, const long to)
{
  assert(to >= 0 && (delay)from < kernel().connection_manager.get_min_delay());
  assert(from < to);

  for (long lag = from; lag < to; ++lag)
  {
    if (S_.r_ == 0)
    {
      // neuron not refractory
      //S_.y3_ = V_.P30_*(S_.y0_ + P_.I_e_) + V_.P33_*S_.y3_ + B_.spikes_.get_value(lag);
      S_.y3_ += S_.y0_ + P_.I_e_ + B_.spikes_.get_value(lag);

      // if we have accumulated spikes from refractory period,
      // add and reset accumulator
      if (P_.with_refr_input_ && S_.refr_spikes_buffer_ != 0.0)
      {
        S_.y3_ += S_.refr_spikes_buffer_;
        S_.refr_spikes_buffer_ = 0.0;
      }

      // lower bound of membrane potential
      S_.y3_ = (S_.y3_ < P_.V_min_ ? P_.V_min_ : S_.y3_);
    }
    else // neuron is absolute refractory
    {
      // read spikes from buffer and accumulate them, discounting
      // for decay until end of refractory period
      if (P_.with_refr_input_)
        S_.refr_spikes_buffer_ += B_.spikes_.get_value(lag) * std::exp(-S_.r_ * Time::get_resolution().get_ms() / P_.tau_m_);
      else
        B_.spikes_.get_value(lag); // clear buffer entry, ignore spike

      --S_.r_;
    }

    // threshold crossing
    if (S_.y3_ >= P_.V_th_)
    {
      S_.r_ = V_.RefractoryCounts_;
      if (P_.subtraction_)
        S_.y3_ = S_.y3_ - P_.V_th_;
      else
        S_.y3_ = P_.V_reset_;

      // EX: must compute spike time
      set_spiketime(Time::step(origin.get_steps() + lag + 1));

      SpikeEvent se;
      kernel().event_delivery_manager.send(*this, se, lag);
    }

    // set new input current
    S_.y0_ = B_.currents_.get_value(lag);

    // voltage logging
    B_.logger_.record_data(origin.get_steps() + lag);
  }
}

void nest::iaf_delta_noleak::handle(SpikeEvent &e)
{
  assert(e.get_delay_steps() > 0);

  // EX: We must compute the arrival time of the incoming spike
  //     explicity, since it depends on delay and offset within
  //     the update cycle.  The way it is done here works, but
  //     is clumsy and should be improved.
  B_.spikes_.add_value(e.get_rel_delivery_steps(kernel().simulation_manager.get_slice_origin()),
                       e.get_weight() * e.get_multiplicity());
}

void nest::iaf_delta_noleak::handle(CurrentEvent &e)
{
  assert(e.get_delay_steps() > 0);

  const double c = e.get_current();
  const double w = e.get_weight();

  // add weighted current; HEP 2002-10-04
  B_.currents_.add_value(e.get_rel_delivery_steps(kernel().simulation_manager.get_slice_origin()),
                         w * c);
}

void nest::iaf_delta_noleak::handle(DataLoggingRequest &e)
{
  B_.logger_.handle(e);
}

} // namespace nest
